Introduction
Today is the first recitation for Module
4 where we put together a lot of the material we’ve learned in the
first 3 modules of this course. Today’s material is on conducting
principal components analysis (PCA) using R, and visualizing the results
with some tools we’ve already learned to use, and some new wrangling and
viz tips along the way.
library(tidyverse) # everything
library(readxl) # reading in excel sheets
library(factoextra) # easy PCA plotting
library(glue) # easy pasting
library(ggrepel) # repelling labels away from their points
library(patchwork) # for combining and arranging plots
Read in data
We will be using data about pizza, which includes data collected
about the nutritional information of 300 different grocery store pizzas,
from 10 brands.
pizza <- read_csv(file = "https://raw.githubusercontent.com/f-imp/Principal-Component-Analysis-PCA-over-3-datasets/master/datasets/Pizza.csv")
## Rows: 300 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): brand
## dbl (8): id, mois, prot, fat, ash, sodium, carb, cal
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
How different are each of the different brands of pizzas analyzed
overall?
1. Run a PCA
Let’s take a look at this new dataset
knitr::kable(head(pizza))
| A |
14069 |
27.82 |
21.43 |
44.87 |
5.11 |
1.77 |
0.77 |
4.93 |
| A |
14053 |
28.49 |
21.26 |
43.89 |
5.34 |
1.79 |
1.02 |
4.84 |
| A |
14025 |
28.35 |
19.99 |
45.78 |
5.08 |
1.63 |
0.80 |
4.95 |
| A |
14016 |
30.55 |
20.15 |
43.13 |
4.79 |
1.61 |
1.38 |
4.74 |
| A |
14005 |
30.49 |
21.28 |
41.65 |
4.82 |
1.64 |
1.76 |
4.67 |
| A |
14075 |
31.14 |
20.23 |
42.31 |
4.92 |
1.65 |
1.40 |
4.67 |
pizza_pca <- prcomp(pizza[,-c(1:2)],
center = TRUE,
scale = TRUE)
We can also look at the output of our PCA in a different way using
the function summary().
summary(pizza_pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.042 1.5134 0.64387 0.3085 0.16636 0.01837 0.003085
## Proportion of Variance 0.596 0.3272 0.05922 0.0136 0.00395 0.00005 0.000000
## Cumulative Proportion 0.596 0.9232 0.98240 0.9960 0.99995 1.00000 1.000000
We can convert this summary into something later usable by extraction
the element importance from
summary(alkaloids_pca) and converting it to a
dataframe.
importance <- summary(pizza_pca)$importance %>%
as.data.frame()
knitr::kable(head(importance))
| Standard deviation |
2.042494 |
1.513426 |
0.6438652 |
0.3085032 |
0.1663641 |
0.0183741 |
0.0030853 |
| Proportion of Variance |
0.595970 |
0.327210 |
0.0592200 |
0.0136000 |
0.0039500 |
0.0000500 |
0.0000000 |
| Cumulative Proportion |
0.595970 |
0.923180 |
0.9824000 |
0.9960000 |
0.9999500 |
1.0000000 |
1.0000000 |
By looking at the summary we can see, for example, that the first two
PCs explain 92.32% of variance.
2. Make a scree plot of the percent variance explained by each
component
Using fviz_eig()
We can do this quickly using fviz_eig().
fviz_eig(pizza_pca)

Manually
If you wanted to make a scree plot manually, you could by plotting
using a wrangled version of the importance dataframe we
made earlier.
# pivot longer
importance_tidy <- importance %>%
rownames_to_column(var = "measure") %>%
pivot_longer(cols = PC1:PC7,
names_to = "PC",
values_to = "value")
# plot
importance_tidy %>%
filter(measure == "Proportion of Variance") %>%
ggplot(aes(x = PC, y = value)) +
geom_col(alpha = 0.1, color = "black") +
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
labs(x = "Principal component",
y = "Percent variance explained",
title = "Pizza scree plot")

3. Make a scores plot of samples, coloring each sample by its
brand
Using fviz_pca_ind()
We can also look at a scores plot using fviz_pca_ind()
where ind means individuals. Here, each point is a sample.
fviz_pca_ind(pizza_pca)

Manually
We want to plot the scores, which are in provided in
pizza_pca$x.
We can convert the list into a dataframe of scores values by using
as.data.frame(). Then we can bind back our relevant
metadata so they’re all together. Note, to use bind_cols()
both datasets need to be in the same order. In this case they are so we
are good.
# create a df of alkaloids_pca$x
scores_raw <- as.data.frame(pizza_pca$x)
# bind meta-data
scores <- bind_cols(pizza[,1], # first columns
scores_raw)
# create objects indicating percent variance explained by PC1 and PC2
PC1_percent <- round((importance[2,1])*100, # index 2nd row, 1st column, times 100
1) # round to 1 decimal
PC2_percent <- round((importance[2,2])*100, 1)
# plot
(scores_plot <- scores %>%
ggplot(aes(x = PC1, y = PC2, fill = brand)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_point(shape = 21, color = "black") +
scale_fill_brewer(palette = "Set3") +
theme_minimal() +
labs(x = glue("PC1: {PC1_percent}%"),
y = glue("PC2: {PC2_percent}%"),
title = "PCA Scores Plot of Pizza Proximate Analysis Across 10 Brands",
fill = "Brand"))

4. Make a loadings plot of samples
Using fviz_pca_var()
We can also look at a loadings plot using fviz_pca_var()
where var means variables. Here, each point is a variable.
fviz_pca_var(pizza_pca)

Manually
We can also make a more customized loadings plot manually using
ggplot and using the dataframe alkaloids_pca$rotation.
# grab raw loadings, without any metadata
loadings_raw <- as.data.frame(pizza_pca$rotation)
# convert rownames to column
loadings <- loadings_raw %>%
rownames_to_column(var = "analysis_type")
# create vector of labels as we want them to appear
analysis_type_labels <- c("Moisture",
"Protein",
"Fat",
"Ash",
"Sodium",
"Carbohydrates",
"Calories")
We can then plot with ggplot like normal.
library(emojifont)
pizza_emoji <- emoji(search_emoji('pizza'))
(loadings_plot <- loadings %>%
ggplot(aes(x = PC1, y = PC2, label = analysis_type_labels)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_point() +
geom_label_repel() +
scale_fill_brewer() +
theme_minimal() +
labs(x = glue("PC1: {PC1_percent}%"),
y = glue("PC2: {PC2_percent}%"),
title = "PCA Loadings Plot for Pizza"))

5. Create either a biplot, or a visualization that shows both your
scores and loadings plot together.
# setting the range of the plot
(scores_plot_ranged <- scores_plot +
coord_cartesian(xlim = c(-3, 5.5), ylim = c(-3, 3.5)))

# what is the ratio of the space on each side of the axis for the scores plot?
(x_ratio <- 3/(3+5.5))
## [1] 0.3529412
(y_ratio <- 3/(3+3.5))
## [1] 0.4615385
# check the ending range for the loadings plot
# 0.7 units looks good for x
x_scalar <- 0.7
# 0.6 units looks good for y
y_scalar <- 0.6
# what should the low range value be so that both plots are equally scaled?
# making the loadings plot match this range
(loadings_plot_ranged <- loadings_plot +
coord_cartesian(xlim = c(-0.3818, x_scalar), ylim = c(-0.5140, y_scalar)))
Plot
scores_plot_ranged + loadings_plot_ranged

Biplot
Using fviz_pca().
You can make a biplot quickly with fviz_pca().
Note, fviz_pca_biplot() and fviz_pca() are the
same.
fviz_pca(pizza_pca)

Instead of making this plot manually, let’s go through how to alter
the existing plot made with fviz_pca(). We can do this
because factoextra creates ggplot objects. To start off, we
need to be using a dataframe that includes our metadata.
# save as a new df
pizza_pca_labelled <- pizza_pca
# assign alkaloid_labels to rownames
rownames(pizza_pca_labelled$rotation) <- analysis_type_labels
# plot
fviz_pca(pizza_pca_labelled, # pca object
label = "var",
repel = TRUE,
geom.var = c("text", "point"),
col.var = "black") +
geom_point(aes(fill = pizza$brand), shape = 21) +
scale_fill_brewer(palette = "Set3") +
theme_minimal() +
labs(x = glue("PC1: {PC1_percent}%"),
y = glue("PC2: {PC2_percent}%"),
title = "PCA Pizza Biplot",
fill = "Brand")

LS0tCnRpdGxlOiAiUHJpbmNpcGFsIENvbXBvbmVudHMgQW5hbHlzaXMg8J+NlSBSZWNpdGF0aW9uIFNvbHV0aW9ucyIKYXV0aG9yOiAiSmVzc2ljYSBDb29wZXJzdG9uZSIKZGF0ZTogIjEwLzI1LzIwMjIiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDUKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgdGhlbWU6IGZsYXRseQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUpCmBgYAoKIyMgSW50cm9kdWN0aW9uCgpUb2RheSBpcyB0aGUgZmlyc3QgcmVjaXRhdGlvbiBmb3IgW01vZHVsZSA0XShfc2l0ZS9tb2R1bGU0Lmh0bWwpIHdoZXJlIHdlIHB1dCB0b2dldGhlciBhIGxvdCBvZiB0aGUgbWF0ZXJpYWwgd2UndmUgbGVhcm5lZCBpbiB0aGUgZmlyc3QgMyBtb2R1bGVzIG9mIHRoaXMgY291cnNlLiBUb2RheSdzIG1hdGVyaWFsIGlzIG9uIGNvbmR1Y3RpbmcgcHJpbmNpcGFsIGNvbXBvbmVudHMgYW5hbHlzaXMgKFBDQSkgdXNpbmcgUiwgYW5kIHZpc3VhbGl6aW5nIHRoZSByZXN1bHRzIHdpdGggc29tZSB0b29scyB3ZSd2ZSBhbHJlYWR5IGxlYXJuZWQgdG8gdXNlLCBhbmQgc29tZSBuZXcgd3JhbmdsaW5nIGFuZCB2aXogdGlwcyBhbG9uZyB0aGUgd2F5LgoKYGBge3IsIGZpZy5hbHQgPSAiQSBwaWN0dXJlIG9mIE5ZIHN0eWxlIGdvb2V5IHBpenphIiwgZmlnLmNhcCA9ICJbU291cmNlXShodHRwczovL3d3dy5hc2VuenlhLmNvbS9ibG9nLzIwMjAvMDMvMjMvbmV3LXlvcmstc3R5bGUtcGl6emEvKSIsIG91dC53aWR0aCA9ICI3MCUiLCBmaWcuYWxpZ24gPSAiY2VudGVyIiwgZWNobyA9IEZBTFNFfQprbml0cjo6aW5jbHVkZV9ncmFwaGljcygiaW1nL3BpenphLnBuZyIpCmBgYAoKYGBge3IgbGlicmFyaWVzLCBtZXNzYWdlID0gRkFMU0UsIHdhcm5pbmcgPSBGQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpICMgZXZlcnl0aGluZwpsaWJyYXJ5KHJlYWR4bCkgIyByZWFkaW5nIGluIGV4Y2VsIHNoZWV0cwpsaWJyYXJ5KGZhY3RvZXh0cmEpICMgZWFzeSBQQ0EgcGxvdHRpbmcKbGlicmFyeShnbHVlKSAjIGVhc3kgcGFzdGluZwpsaWJyYXJ5KGdncmVwZWwpICMgcmVwZWxsaW5nIGxhYmVscyBhd2F5IGZyb20gdGhlaXIgcG9pbnRzCmxpYnJhcnkocGF0Y2h3b3JrKSAjIGZvciBjb21iaW5pbmcgYW5kIGFycmFuZ2luZyBwbG90cwpgYGAKCiMjIyMgUmVhZCBpbiBkYXRhCgpXZSB3aWxsIGJlIHVzaW5nIGRhdGEgYWJvdXQgcGl6emEsIHdoaWNoIGluY2x1ZGVzIGRhdGEgY29sbGVjdGVkIGFib3V0IHRoZSBudXRyaXRpb25hbCBpbmZvcm1hdGlvbiBvZiAzMDAgZGlmZmVyZW50IGdyb2Nlcnkgc3RvcmUgcGl6emFzLCBmcm9tIDEwIGJyYW5kcy4KCmBgYHtyfQpwaXp6YSA8LSByZWFkX2NzdihmaWxlID0gImh0dHBzOi8vcmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbS9mLWltcC9QcmluY2lwYWwtQ29tcG9uZW50LUFuYWx5c2lzLVBDQS1vdmVyLTMtZGF0YXNldHMvbWFzdGVyL2RhdGFzZXRzL1BpenphLmNzdiIpCmBgYAoKPiBIb3cgZGlmZmVyZW50IGFyZSBlYWNoIG9mIHRoZSBkaWZmZXJlbnQgYnJhbmRzIG9mIHBpenphcyBhbmFseXplZCBvdmVyYWxsPyAKCiMjIDEuIFJ1biBhIFBDQQoKTGV0J3MgdGFrZSBhIGxvb2sgYXQgdGhpcyBuZXcgZGF0YXNldApgYGB7ciBoZWFkfQprbml0cjo6a2FibGUoaGVhZChwaXp6YSkpCmBgYAoKYGBge3IgcHJjb21wfQpwaXp6YV9wY2EgPC0gcHJjb21wKHBpenphWywtYygxOjIpXSwKICAgICAgICAgICAgICAgICAgICBjZW50ZXIgPSBUUlVFLAogICAgICAgICAgICAgICAgICAgIHNjYWxlID0gVFJVRSkKYGBgCgpXZSBjYW4gYWxzbyBsb29rIGF0IHRoZSBvdXRwdXQgb2Ygb3VyIFBDQSBpbiBhIGRpZmZlcmVudCB3YXkgdXNpbmcgdGhlIGZ1bmN0aW9uIGBzdW1tYXJ5KClgLgpgYGB7ciBzdW1tYXJ5IHBjYX0Kc3VtbWFyeShwaXp6YV9wY2EpIApgYGAKCldlIGNhbiBjb252ZXJ0IHRoaXMgc3VtbWFyeSBpbnRvIHNvbWV0aGluZyBsYXRlciB1c2FibGUgYnkgZXh0cmFjdGlvbiB0aGUgZWxlbWVudCBgaW1wb3J0YW5jZWAgZnJvbSBgc3VtbWFyeShhbGthbG9pZHNfcGNhKWAgYW5kIGNvbnZlcnRpbmcgaXQgdG8gYSBkYXRhZnJhbWUuCgpgYGB7ciBleHRyYWN0IGltcG9ydGFuY2V9CmltcG9ydGFuY2UgPC0gc3VtbWFyeShwaXp6YV9wY2EpJGltcG9ydGFuY2UgJT4lCiAgYXMuZGF0YS5mcmFtZSgpCgprbml0cjo6a2FibGUoaGVhZChpbXBvcnRhbmNlKSkKYGBgCgpCeSBsb29raW5nIGF0IHRoZSBzdW1tYXJ5IHdlIGNhbiBzZWUsIGZvciBleGFtcGxlLCB0aGF0IHRoZSBmaXJzdCB0d28gUENzIGV4cGxhaW4gYHIgcm91bmQoKGltcG9ydGFuY2VbMywyXSkqMTAwLCAyKWAlIG9mIHZhcmlhbmNlLgoKIyMgMi4gTWFrZSBhIHNjcmVlIHBsb3Qgb2YgdGhlIHBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkIGJ5IGVhY2ggY29tcG9uZW50CgojIyMgVXNpbmcgYGZ2aXpfZWlnKClgCldlIGNhbiBkbyB0aGlzIHF1aWNrbHkgdXNpbmcgW2Bmdml6X2VpZygpYF0oaHR0cHM6Ly9ycGtncy5kYXRhbm92aWEuY29tL2ZhY3RvZXh0cmEvcmVmZXJlbmNlL2VpZ2VudmFsdWUuaHRtbCkuCmBgYHtyIGZ2aXpfZWlnfQpmdml6X2VpZyhwaXp6YV9wY2EpCmBgYAoKIyMjIE1hbnVhbGx5CklmIHlvdSB3YW50ZWQgdG8gbWFrZSBhIHNjcmVlIHBsb3QgbWFudWFsbHksIHlvdSBjb3VsZCBieSBwbG90dGluZyB1c2luZyBhIHdyYW5nbGVkIHZlcnNpb24gb2YgdGhlIGBpbXBvcnRhbmNlYCBkYXRhZnJhbWUgd2UgbWFkZSBlYXJsaWVyLgoKYGBge3IgaW1wb3J0YW5jZSB3cmFuZ2xpbmcgbWFudWFsIHBsb3QgMX0KIyBwaXZvdCBsb25nZXIKaW1wb3J0YW5jZV90aWR5IDwtIGltcG9ydGFuY2UgJT4lCiAgcm93bmFtZXNfdG9fY29sdW1uKHZhciA9ICJtZWFzdXJlIikgJT4lCiAgcGl2b3RfbG9uZ2VyKGNvbHMgPSBQQzE6UEM3LAogICAgICAgICAgICAgICBuYW1lc190byA9ICJQQyIsCiAgICAgICAgICAgICAgIHZhbHVlc190byA9ICJ2YWx1ZSIpCgojIHBsb3QKaW1wb3J0YW5jZV90aWR5ICU+JQogIGZpbHRlcihtZWFzdXJlID09ICJQcm9wb3J0aW9uIG9mIFZhcmlhbmNlIikgJT4lCiAgZ2dwbG90KGFlcyh4ID0gUEMsIHkgID0gdmFsdWUpKSArCiAgZ2VvbV9jb2woYWxwaGEgPSAwLjEsIGNvbG9yID0gImJsYWNrIikgKwogIHNjYWxlX3lfY29udGludW91cyhsYWJlbHMgPSBzY2FsZXM6OnBlcmNlbnQpICsKICB0aGVtZV9taW5pbWFsKCkgKwogIGxhYnMoeCA9ICJQcmluY2lwYWwgY29tcG9uZW50IiwKICAgICAgIHkgPSAiUGVyY2VudCB2YXJpYW5jZSBleHBsYWluZWQiLAogICAgICAgdGl0bGUgPSAiUGl6emEgc2NyZWUgcGxvdCIpCmBgYAoKIyMgMy4gTWFrZSBhIHNjb3JlcyBwbG90IG9mIHNhbXBsZXMsIGNvbG9yaW5nIGVhY2ggc2FtcGxlIGJ5IGl0cyBicmFuZAoKIyMjIFVzaW5nIGBmdml6X3BjYV9pbmQoKWAKCldlIGNhbiBhbHNvIGxvb2sgYXQgYSBzY29yZXMgcGxvdCB1c2luZyBbYGZ2aXpfcGNhX2luZCgpYF0oaHR0cHM6Ly9ycGtncy5kYXRhbm92aWEuY29tL2ZhY3RvZXh0cmEvcmVmZXJlbmNlL2Z2aXpfcGNhLmh0bWwpIHdoZXJlIGluZCBtZWFucyBpbmRpdmlkdWFscy4gSGVyZSwgZWFjaCBwb2ludCBpcyBhIHNhbXBsZS4KYGBge3Igc2NvcmVzIGZ2aXpfcGNhX2luZH0KZnZpel9wY2FfaW5kKHBpenphX3BjYSkKYGBgCgojIyMjIE1hbnVhbGx5CgpXZSB3YW50IHRvIHBsb3QgdGhlIHNjb3Jlcywgd2hpY2ggYXJlIGluIHByb3ZpZGVkIGluIGBwaXp6YV9wY2EkeGAuCgpXZSBjYW4gY29udmVydCB0aGUgbGlzdCBpbnRvIGEgZGF0YWZyYW1lIG9mIHNjb3JlcyB2YWx1ZXMgYnkgdXNpbmcgYGFzLmRhdGEuZnJhbWUoKWAuIFRoZW4gd2UgY2FuIGJpbmQgYmFjayBvdXIgcmVsZXZhbnQgbWV0YWRhdGEgc28gdGhleSdyZSBhbGwgdG9nZXRoZXIuIE5vdGUsIHRvIHVzZSBgYmluZF9jb2xzKClgIGJvdGggZGF0YXNldHMgbmVlZCB0byBiZSBpbiB0aGUgc2FtZSBvcmRlci4gSW4gdGhpcyBjYXNlIHRoZXkgYXJlIHNvIHdlIGFyZSBnb29kLgpgYGB7ciBzY29yZXMgd3JhbmdsaW5nfQojIGNyZWF0ZSBhIGRmIG9mIGFsa2Fsb2lkc19wY2EkeApzY29yZXNfcmF3IDwtIGFzLmRhdGEuZnJhbWUocGl6emFfcGNhJHgpCgojIGJpbmQgbWV0YS1kYXRhCnNjb3JlcyA8LSBiaW5kX2NvbHMocGl6emFbLDFdLCAjIGZpcnN0IGNvbHVtbnMKICAgICAgICAgICAgICAgICAgICBzY29yZXNfcmF3KQpgYGAKCgpgYGB7ciBzY29yZXMgcGxvdHRpbmcgMn0KIyBjcmVhdGUgb2JqZWN0cyBpbmRpY2F0aW5nIHBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkIGJ5IFBDMSBhbmQgUEMyClBDMV9wZXJjZW50IDwtIHJvdW5kKChpbXBvcnRhbmNlWzIsMV0pKjEwMCwgIyBpbmRleCAybmQgcm93LCAxc3QgY29sdW1uLCB0aW1lcyAxMDAKICAgICAgICAgICAgICAgICAgICAgMSkgIyByb3VuZCB0byAxIGRlY2ltYWwKUEMyX3BlcmNlbnQgPC0gcm91bmQoKGltcG9ydGFuY2VbMiwyXSkqMTAwLCAxKSAKCiMgcGxvdAooc2NvcmVzX3Bsb3QgPC0gc2NvcmVzICU+JQogIGdncGxvdChhZXMoeCA9IFBDMSwgeSA9IFBDMiwgZmlsbCA9IGJyYW5kKSkgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiKSArCiAgZ2VvbV9wb2ludChzaGFwZSA9IDIxLCBjb2xvciA9ICJibGFjayIpICsKICBzY2FsZV9maWxsX2JyZXdlcihwYWxldHRlID0gIlNldDMiKSArCiAgdGhlbWVfbWluaW1hbCgpICsKICBsYWJzKHggPSBnbHVlKCJQQzE6IHtQQzFfcGVyY2VudH0lIiksIAogICAgICAgeSA9IGdsdWUoIlBDMjoge1BDMl9wZXJjZW50fSUiKSwgCiAgICAgICB0aXRsZSA9ICJQQ0EgU2NvcmVzIFBsb3Qgb2YgUGl6emEgUHJveGltYXRlIEFuYWx5c2lzIEFjcm9zcyAxMCBCcmFuZHMiLAogICAgICAgZmlsbCA9ICJCcmFuZCIpKQpgYGAKCiMjIDQuIE1ha2UgYSBsb2FkaW5ncyBwbG90IG9mIHNhbXBsZXMKCiMjIyBVc2luZyBgZnZpel9wY2FfdmFyKClgCgpXZSBjYW4gYWxzbyBsb29rIGF0IGEgbG9hZGluZ3MgcGxvdCB1c2luZyBbYGZ2aXpfcGNhX3ZhcigpYF0oaHR0cHM6Ly9ycGtncy5kYXRhbm92aWEuY29tL2ZhY3RvZXh0cmEvcmVmZXJlbmNlL2Z2aXpfcGNhLmh0bWwpIHdoZXJlIHZhciBtZWFucyB2YXJpYWJsZXMuIEhlcmUsIGVhY2ggcG9pbnQgaXMgYSB2YXJpYWJsZS4KYGBge3IgbG9hZGluZ3MgZnZpel9wY2FfdmFyfQpmdml6X3BjYV92YXIocGl6emFfcGNhKQpgYGAKCiMjIyBNYW51YWxseQoKV2UgY2FuIGFsc28gbWFrZSBhIG1vcmUgY3VzdG9taXplZCBsb2FkaW5ncyBwbG90IG1hbnVhbGx5IHVzaW5nIGdncGxvdCBhbmQgdXNpbmcgdGhlIGRhdGFmcmFtZSBgYWxrYWxvaWRzX3BjYSRyb3RhdGlvbmAuCgpgYGB7ciBsb2FkaW5ncyB3cmFuZ2xpbmd9CiMgZ3JhYiByYXcgbG9hZGluZ3MsIHdpdGhvdXQgYW55IG1ldGFkYXRhCmxvYWRpbmdzX3JhdyA8LSBhcy5kYXRhLmZyYW1lKHBpenphX3BjYSRyb3RhdGlvbikKCiMgY29udmVydCByb3duYW1lcyB0byBjb2x1bW4KbG9hZGluZ3MgPC0gbG9hZGluZ3NfcmF3ICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbih2YXIgPSAiYW5hbHlzaXNfdHlwZSIpCgojIGNyZWF0ZSB2ZWN0b3Igb2YgbGFiZWxzIGFzIHdlIHdhbnQgdGhlbSB0byBhcHBlYXIKYW5hbHlzaXNfdHlwZV9sYWJlbHMgPC0gYygiTW9pc3R1cmUiLAogICAgICAgICAgICAgICAgICAgICAgICAgICJQcm90ZWluIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAiRmF0IiwKICAgICAgICAgICAgICAgICAgICAgICAgICAiQXNoIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAiU29kaXVtIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAiQ2FyYm9oeWRyYXRlcyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIkNhbG9yaWVzIikKYGBgCgpXZSBjYW4gdGhlbiBwbG90IHdpdGggZ2dwbG90IGxpa2Ugbm9ybWFsLgoKYGBge3IgbG9hZGluZ3N9CmxpYnJhcnkoZW1vamlmb250KQpwaXp6YV9lbW9qaSA8LSBlbW9qaShzZWFyY2hfZW1vamkoJ3BpenphJykpCgoobG9hZGluZ3NfcGxvdCA8LSBsb2FkaW5ncyAlPiUKICBnZ3Bsb3QoYWVzKHggPSBQQzEsIHkgPSBQQzIsIGxhYmVsID0gYW5hbHlzaXNfdHlwZV9sYWJlbHMpKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCwgbGluZXR5cGUgPSAiZGFzaGVkIikgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGFiZWxfcmVwZWwoKSArCiAgc2NhbGVfZmlsbF9icmV3ZXIoKSArCiAgdGhlbWVfbWluaW1hbCgpICsKICBsYWJzKHggPSBnbHVlKCJQQzE6IHtQQzFfcGVyY2VudH0lIiksIAogICAgICAgeSA9IGdsdWUoIlBDMjoge1BDMl9wZXJjZW50fSUiKSwgCiAgICAgICB0aXRsZSA9ICJQQ0EgTG9hZGluZ3MgUGxvdCBmb3IgUGl6emEiKSkKYGBgCgojIyA1LiBDcmVhdGUgZWl0aGVyIGEgYmlwbG90LCBvciBhIHZpc3VhbGl6YXRpb24gdGhhdCBzaG93cyBib3RoIHlvdXIgc2NvcmVzIGFuZCBsb2FkaW5ncyBwbG90IHRvZ2V0aGVyLgoKYGBge3IgcGF0Y2h3b3JrIHJlLXJhbmdpbmcgYW5kIHBsb3R9CiMgc2V0dGluZyB0aGUgcmFuZ2Ugb2YgdGhlIHBsb3QKKHNjb3Jlc19wbG90X3JhbmdlZCA8LSBzY29yZXNfcGxvdCArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC0zLCA1LjUpLCB5bGltID0gYygtMywgMy41KSkpCgojIHdoYXQgaXMgdGhlIHJhdGlvIG9mIHRoZSBzcGFjZSBvbiBlYWNoIHNpZGUgb2YgdGhlIGF4aXMgZm9yIHRoZSBzY29yZXMgcGxvdD8KKHhfcmF0aW8gPC0gMy8oMys1LjUpKQooeV9yYXRpbyA8LSAzLygzKzMuNSkpCgojIGNoZWNrIHRoZSBlbmRpbmcgcmFuZ2UgZm9yIHRoZSBsb2FkaW5ncyBwbG90CiMgMC43IHVuaXRzIGxvb2tzIGdvb2QgZm9yIHgKeF9zY2FsYXIgPC0gMC43CiMgMC42IHVuaXRzIGxvb2tzIGdvb2QgZm9yIHkKeV9zY2FsYXIgPC0gMC42CgojIHdoYXQgc2hvdWxkIHRoZSBsb3cgcmFuZ2UgdmFsdWUgYmUgc28gdGhhdCBib3RoIHBsb3RzIGFyZSBlcXVhbGx5IHNjYWxlZD8KCiMgbWFraW5nIHRoZSBsb2FkaW5ncyBwbG90IG1hdGNoIHRoaXMgcmFuZ2UKKGxvYWRpbmdzX3Bsb3RfcmFuZ2VkIDwtIGxvYWRpbmdzX3Bsb3QgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtMC4zODE4LCB4X3NjYWxhciksIHlsaW0gPSBjKC0wLjUxNDAsIHlfc2NhbGFyKSkpCmBgYApQbG90CmBgYHtyIHBhdGNod29yayByZXJhbmdlZCwgZmlnLndpZHRoID0gMTB9CnNjb3Jlc19wbG90X3JhbmdlZCArIGxvYWRpbmdzX3Bsb3RfcmFuZ2VkCmBgYAoKIyMgQmlwbG90CgojIyMgVXNpbmcgYGZ2aXpfcGNhKClgLgoKWW91IGNhbiBtYWtlIGEgYmlwbG90IHF1aWNrbHkgd2l0aCBbYGZ2aXpfcGNhKClgXShodHRwczovL3Jwa2dzLmRhdGFub3ZpYS5jb20vZmFjdG9leHRyYS9yZWZlcmVuY2UvZnZpel9wY2EuaHRtbCkuIE5vdGUsIGBmdml6X3BjYV9iaXBsb3QoKWAgYW5kIGBmdml6X3BjYSgpYCBhcmUgdGhlIHNhbWUuCmBgYHtyIGZ2aXpfcGNhfQpmdml6X3BjYShwaXp6YV9wY2EpCmBgYAoKSW5zdGVhZCBvZiBtYWtpbmcgdGhpcyBwbG90IG1hbnVhbGx5LCBsZXQncyBnbyB0aHJvdWdoIGhvdyB0byBhbHRlciB0aGUgZXhpc3RpbmcgcGxvdCBtYWRlIHdpdGggYGZ2aXpfcGNhKClgLiBXZSBjYW4gZG8gdGhpcyBiZWNhdXNlIGBmYWN0b2V4dHJhYCBjcmVhdGVzIGdncGxvdCBvYmplY3RzLiBUbyBzdGFydCBvZmYsIHdlIG5lZWQgdG8gYmUgdXNpbmcgYSBkYXRhZnJhbWUgdGhhdCBpbmNsdWRlcyBvdXIgbWV0YWRhdGEuCgpgYGB7ciBmdml6X3BjYSByZWxhYmVsbGVkIHBsb3QgMn0KIyBzYXZlIGFzIGEgbmV3IGRmCnBpenphX3BjYV9sYWJlbGxlZCA8LSBwaXp6YV9wY2EKCiMgYXNzaWduIGFsa2Fsb2lkX2xhYmVscyB0byByb3duYW1lcwpyb3duYW1lcyhwaXp6YV9wY2FfbGFiZWxsZWQkcm90YXRpb24pIDwtIGFuYWx5c2lzX3R5cGVfbGFiZWxzCgojIHBsb3QKZnZpel9wY2EocGl6emFfcGNhX2xhYmVsbGVkLCAjIHBjYSBvYmplY3QKICAgICAgICAgbGFiZWwgPSAidmFyIiwKICAgICAgICAgcmVwZWwgPSBUUlVFLAogICAgICAgICBnZW9tLnZhciA9IGMoInRleHQiLCAicG9pbnQiKSwKICAgICAgICAgY29sLnZhciA9ICJibGFjayIpICsKICBnZW9tX3BvaW50KGFlcyhmaWxsID0gcGl6emEkYnJhbmQpLCBzaGFwZSA9IDIxKSArCiAgc2NhbGVfZmlsbF9icmV3ZXIocGFsZXR0ZSA9ICJTZXQzIikgKwogIHRoZW1lX21pbmltYWwoKSArCiAgbGFicyh4ID0gZ2x1ZSgiUEMxOiB7UEMxX3BlcmNlbnR9JSIpLCAKICAgICAgIHkgPSBnbHVlKCJQQzI6IHtQQzJfcGVyY2VudH0lIiksIAogICAgICAgdGl0bGUgPSAiUENBIFBpenphIEJpcGxvdCIsCiAgICAgICBmaWxsID0gIkJyYW5kIikKYGBgCgoK